arXiv:1508.00563vl [astro-ph.EP] 3 Aug 2015 


Draft version August 5, 2015 

Preprint typeset using DTeX style emulateapj v. 5/2/11 


IMPROVED ALGORITHMS FOR RADAR-BASED RECONSTRUCTION OF ASTEROID SHAPES 

Adam H. Greenberg and Jean-Luc Margot 
University California, Los Angeles 
Draft version August 5, 2015 

ABSTRACT 

We describe our implementation of a global-parameter optimizer and Square Root Information 
Filter (SRIF) into the asteroid-modelling software shape. We compare the performance of our new 
optimizer with that of the existing sequential optimizer when operating on various forms of simulated 
data and actual asteroid radar data. In all cases, the new implementation performs substantially 
better than its predecessor: it converges faster, produces shape models that are more accurate, and 
solves for spin axis orientations more reliably. We discuss potential future changes to improve shape’s 
fitting speed and accuracy. 

Subject headings: asteroids, 2000 ET70, radar, shape, model, optimization, SRIF 


I. INTRODUCTION 

Earth-based radar is a powerful tool for gathering in¬ 
formation about bodies in the Solar System. Radar ob¬ 
servations can dramatically improve the determination 
of the physical properties and orbital elements of small 
bodies (such as asteroids and comets). An important 
development in the past two decades has been the for¬ 
mulation and implementation of algorith ms for asteroid 
shape reconstruction based on radar data (|Hudson|I993 


Hudson fc: Ostro|I994[|Ostro et al. 


udson V Ostro 


I995p. This problem is not trivial because it requires the 

joint estimation of the spin state and shape of the aster¬ 
oid. Because of the nature of radar data, recovery of the 
spin state depends on knowledge of the shape and vice 
versa. Even with perfect spin state information, certain 
peculiarities of radar images (such as the two-to-one or 
several-to-one mapping between surface elements on the 
object and pixels within the radar image) mak e recovery 
of the physical shape challenging (Ostro||I993). This is a 
computationally intensive problem, potentially involving 
hundreds to thousands of free parameters and millions of 
data points. 

Despite the computational cost, astronomers are keen 
on deriving shape and spin information from asteroid 
radar images. The most compelling reason to do so is 
the fact that radar is the only Earth-based technique 
that can produce detailed three-dimensional information 
of near-Earth objects. This is possible because radar in¬ 
struments achieve spatial resolutions that dramatically 
surpass the diffraction limit. In other words, radar in¬ 
struments can resolve objects substantially smaller than 
the beamwidth of the antenna used to obtain the images. 
For example, the Arecibo telescope, the primary instru¬ 
ment used for the data presented in this paper, has a 
beamwidth of ^2 arcminutes at the nominal 2380 MHz 
frequency of the radar. Yet observers can easily gather 
shape information to an accuracy of decameters for ob¬ 
jects several millions of kilometers from Earth, achieving 
an effective spatial resolution of ^I milliarcsecond. 

Radar has other advantages as well. Unlike most ob¬ 
servational techniques inside the Solar System, radar 
does not rely on any external sources of light, be it re¬ 
flected sunlight, transmitted starlight, or thermal emis¬ 


sion. This human-controlled illumination allows for 
greater flexibility with respect to the observations. In 
addition, because of the wavelengths involved, radar ob¬ 
servations can be performed during the day, further en¬ 
hancing this flexibility. Radar also has the ability to 
probe an object’s sub-surface properties, which can give 
important information about the object such as poros¬ 
ity, surface and sub-surface dielectric constant, and the 
presence of near-surface ice. 

Asteroid shape data are important for various reasons. 
For certain asteroids, reliable determination of an or¬ 
bital future cannot be determined without shape and 
spin information. The Yarkovsky effect, for example, 
can change an asteroid’s semi-maio r axis at a rate of 

10“^ AU/My for km-sized objects jVokrouhlicky et al ' 


2000 Bottke et al. 2006 INugent et al. 2012L This ef¬ 


fect occurs because the rotating body absorbs sunlight 
and then re-emits that light in a non-sunward direction, 
resulting in a gentle perturbation to the asteroid’s or¬ 
bit. The Yarkovsky effect is greatly dependent on the 
shape of the object, since re-emission of absorbed sun¬ 
light is a surface phenomenon. It is responsible for the 
largest source of uncertainty in trajectory predictions for 
near-Earth asteroids (NEAs) with sizes under 2 km, and 
it must be ta ken into account when evaluating impact 
probabilities (|Giorgini et ah 2002 Chesley et al. 2014 


Farnocchia et al. pOTTl r 


Knowledge of the shape also provides clues about the 
formation and interaction history of asteroids. For ex¬ 
ample, radar-derived shapes of asteroids have been in¬ 
strumental in identifying binary asteroids and contact 
binaries, which repre sent ^16% and ^10% of the pop¬ 


ulatio n, respectively (Margot et al. 2002 Benner et al. 


2008). They have also provided strong evidence that 


INF A binaries form by a spin- up and mass sh edding pro¬ 


cess (Margot et al. 2002 Ostro et al. 2006). For single 


asteroids, knowledge of the morphology guides interpre¬ 
tation of the collisional history and surface modification 
processes. 

Shapes also affect spin evolution during two-body in¬ 
teractions (e.g., torques during close planetary encoun¬ 
ters) and orbital evolution of binary NEAs (e.g., tidal, 
gravitational, and non-gravitational in^ t eractions between 
components) (e.g., Margot et al.[|2QQ2 Ostro et ar]|2QQ6 
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Figure 1. A time-series of range-Doppler images of the asteroid 2000 ET70 (|Naidu et al.|2013|), starting in the top left and proceeding to 
the right. The epochs of consecutive images are separated by 18 minutes. Distance trom the observer increases downwards, and Doppler 
increases to the right. 


Scheeres et al.||2QQ6| Cnk fc NesvQrny||2QlQ 

et al. 201^ Naidu &: MargQt||2Q15p . Finally, s 


Jacobson 


napes are 

needed when calculating the gravity environment near 


operations ( 

Eujiwara et al. 

|2^ 

i Naidu et al.|2013 

Nolan 

et al.||2013 

Takahashi & Scheeres||2U14|). 


radar data is equally valuable. In contrast to light curve 
period determinations, which are neither sidereal nor 
synodic, the radar-based measurements yield sidereal pe¬ 
riods. These estimates are needed to test the agree¬ 
ment between physical theories and observations, e.g.. 


lor et al. 112007 iLowry et al.||2007 

) and subsequent shape 

evolution (e.g.,|Harris et aL|: 

2(m 

Eahnestock & Scheeres 


2009). Proper modeling ot the Yarkovsky perturbations 


to an asteroid’s heliocentric orbit or to the evolution 


of binary orbits (e.g., Margot et al. 2015) also require 
knowledge of the spin state, h'inahy, important insights 
can be gained about asteroid physical properties and col- 
lisional evolution from the spin distributions of b oth reg¬ 
ular rotator s and non-principal-axis rotators (e.g., Pravec 
etaLl[2002l ). 


2. CURRENT METHOD 


to the model parameters in an attempt to minimize the 
sum of squares of residuals, shape uses increasing model 
complexity to build up a representation for the asteroid, 
from a basic ellipsoid model to capture gross features, 
to a spherical harmonic model whi ch can represent hner 
surface elements (See section 5.1) and hnahy a model 
based on contiguous triangular facets (hereafter vertex 
model). The spin state is generally estimated in the early 
stages of the shape htting - this is normally done by using 
trial values of the spin state while simultaneously htting 
for the shape itself. 

shape currently uses a Sequential Parameter Fit (SPF) 
mechanism to adjust the model following a comparison 
between the model projection and the radar observables. 
SPF minimizes using a “bracket and Brent” method 
( Press et al.|1992 ) - for each iteration, this process mini- 
mizes for variations in that individual parameter only, 
while ah other parameters are held constant. This pro¬ 
cess is not only slow, but it also does not guarantee con¬ 
vergence on a global minimum, or even the nearest local 
minimum, because minimization always progresses along 
a single parameter axis at a time. We have worked to¬ 
wards replacing the SPF currently implemented in shape 
with a modihed Squa re R oot Information Filter (SRIF), 
as outlined in section [3Tl 


Asteroid shapes and spin states ar e currently modeled 


using the shap e software package ( Hudson] 1993 


Ma- 


gri et al.|2007). shape takes a model for the asteroid, 
which is based on both shape and spin parameters, as 
well as scattering behavior, and projects that model into 
the same space as that of the radar observables. This 
space, called the range-Doppler space, has dimensions of 
range and line-of-sight velocity (hgurel^. shape can also 
handle optical lightcurve data when fitting for asteroid 
shapes, but we did not use this capability in this paper. 

shape then compares the mapping of the model into 
this space to the radar observables, and makes changes 


3. SOLUTION VIA NORMAL EQUATIONS 


Before detailing the mechanics of the SRIF, it is worth 
discussing the Normal Equations Method (NEM), to 


which SRIE is related ( 

Press et al. 

11992). A classical 

NEM minimizes the weighted residua 

Is between a model 


and data with noise assumed to be Gaussian by deter¬ 
mining the direction in parameter space in which is 
decreasing fastest. Specihcahy, suppose one has a set of 
m observables, with weights that are the diagonal ele¬ 
ments of an m X m matrix kF, and a model function /(x), 
where x is an n-dimensional parameter vector. Assum- 







































































Improved Algorithms for Radar-Based Reconstruction of Asteroid Shapes 


3 


ing independent data points with Gaussian-distributed 
errors, the probability of the model matching the data is 
given by 


p{f{x) I i*) (xp{z I /(f)) oc exp{-^R'^WR), 

where R = z — f{x). Therefore maximizing the model 
probability is the same as minimizing the value 

X^{x)=WWR. 

Perturbing f by some amount, 5x, and minimizing x^(f) 
over 5x yields 

(AWA)fe = AWR, 

where 



Thus, changing one’s parameter vector by 

51e = {A^WA)-^A^WR ( 1 ) 

yields a decrease in x^(x). For non-linear systems, this 
process is repeated multiple times until the change in 
from one iteration to the next has passed below a fidu¬ 
cial fraction. Equation is also known as the weighted 
normal equation. 

A major issue with NEM is the computation of the 
inverse of the matrix A^WA. This matrix has v? el¬ 
ements and thus can be quite large for a model with 
many parameters. In addition, numerical stability can 
be a serious issue - A'^WA may be ill-conditioned, and 
thus taking the inverse can result in numerical errors (see 
Appendix). 

One way to quantify the issue of numerical stabil¬ 
ity is by using the condition number where 

ti{M) = ||M|| * ||M“^||. A smaller /^(M) corresponds 
to a better conditioned matrix M, meaning that fewer 
errors will accrue in the calculation of M~^. 

Since 

/^(AWA) (X K{Af 

and 

k{A) > 1 

for non-orthogonal matrix A, the classical NEM increases 
the risk of numerical instabilities 

Einally, for problems involving a very large number of 
observations and model parameters, even the calculation 
of {A'^WA) is non-trivial, as this matrix multiplication 
scales like im?n. The number of observations needed for 
an asteroid shape reconstruction typically number in the 
millions, with potentially 10^ — 10^ free parameters. 


3.1. Square Root Information Filter 
The Square Root Information Eilter ( SRIE) was orig¬ 


i nally developed by Bierman in 1977 (|Bierman, G. J. 
(1977), Lawson & Hanson (1995)). The algorithm mini- 
mizes for time series data with Gaussian errors, and 
is inspired by the Kalman filter algorithm. SRIE is more 
stable, more accurate, and faster than the algorithm cur¬ 
rently used in shape. SRIE is also more numerically 


stable (and, in some cases faster) than a solution via 
normal equations. Our implementation of SRIE includes 
some changes to the original algorithm, which will be 
discussed in section UJ 

SRIE gets around all the problems described above by 
utili zing matrix square root s and Householder operations 
(see Bierman, G. J. (1977), pg. 59) to increase the nu- 
merical stability when determining 6x. Instead of mini¬ 
mizing x^, SRIE minimizes 


Q = {xT- =\\w--R\l 


where IF 2 is the square root of the matrix IF, defined 
such that 

IF = IF^IFF 


In general, the square root of a matrix is multivalued, 
however since IF is positive-semidefinite, all square roots 
are real. We select the positive root by convention. 

Then, along similar lines as NEM, a change of 6x is 
introduced to the parameter vector F, and Q' = Q{x-\-Sx) 
is minimized over this change. 

Q' is smallest when 

\\W^R{xFSx)\\ ^ \\W^R{x) F ASx)\\ 

= ||IF^R(f) + IF^Afe|| 

is minimized. 

A matrix H is defined such that HW 2 A is upper tri¬ 
angular. H is orthogonal and can be generated by the 
product of n Householder operation matrices. Note that 
the orthogonality of H guarantees that 

||IF^R(f) + IF^Afe|| = ||i7(IF^R(f) + IF^Afe)|| 

= \\HW"^R{x) F HW^A5lr\\. 

Erom the definition of i7, HW^A can be rewritten as 


HW^A 



where A' is an n x n upper-triangular matrix, and Z is 
the {m — n) X n zero-matrix. Then, rewriting 

HWiR{x) = ( ) > 

where R'^ and R'^ are m x 1 and {n — m) x 1 arrays, 
respectively, yields 


f R'^ + A'sic \ 



\ F Z5x J 


V 4 ) 


This is clearly minimized over Sx when 
R'^ = -A'Sx 


or 

Sx = -A'-^R'^. (2) 

Since A' is upper triangular, its inverse can be easily 
calculated, and singularity can be trivially detected. Eur- 
thermore, the condition number of the inverted matrix is 
proportional to /^(A), as opposed to i<i{A)‘^ in the NEM 
case. 
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Finally, note that A'^WA is never calculated, which, 
as mentioned in section is a computationally intensive 
calculation. 

4. ADDITIONS TO SRIF 
4.1. Optimizations 

The number of operations necessary to generate the 
Householder matrix H grows as 0(n^(m — n)), where 
the number of observations m always exceeds the num¬ 
ber of parameters n. Although this growth profile is 
favorable with respect to m when compared to that of 
NEM (O(m^n)), it becomes problematic for high reso¬ 
lution models (large n). To maintain good performance 
in large n situations, we have implemented three main 
optimizations to the standard SRIF. 

Our first addition is to run the matrix triangulariza- 
tion simultaneously on multiple cores, which results in 
a significant speed-up. Note that although Householder 
matrices are generated iteratively, any given iteration k 
requires n — k column-wise operations, and each of these 
operations are independent from each-other. Therefore, 
the Householder matrix calculations can be done in a 
thread-safe manner. 

The second addition we made to the standard SRIF 
is the inclusion of a secondary minimization for the 
scaling of Sx, so that 


4.2. Penalty functions 

The SPF routine can currently fit models to data while 
taking into account a suite of “penalty functions” that 
favor models with desirable properties. In a way, these 
penalty functions serve to make the fit operate in a more 
global context - there may be a local minimum in y^- 
space towards which the fitting algorithm would want to 
tend, but that minimum can be ruled out a priori thanks 
to physical considerations. These penalty functions in¬ 
clude limits on ellipsoid axis ratios to avoid absurdly 
elongated or flattened shapes, constraints on shape con¬ 
cavities to avoid unrealistic surface topographies, and 
limits on the model center of mass distance from the im¬ 
age center, to name a few. We have implemented these 
penalty functions in the SRIF framework by redefining 
the residual vector as 

R'> = ( 

\ Pw J 

where 

PlXWi \ 

Pn xwn / 

for which {pi} , {wi} are the set of penalty functions and 
penalty weights, respectively, and 


Q' = Q{x + a6x) 

is minimized over a. This minimization is done with 
an eleven point grid search for a, from a = 10“^ to 
a = 10^’^. The additional minimization adds a trivial 
additional computation cost to the overall minimization 
of y^, but allows for faster convergence, and the possi¬ 
bility of skipping over local minima in the y^-space. 

The final change we made to the underlying SRIF algo¬ 
rithm also granted the largest speed improvement. Even 
with the optimizations described above, the 0{n^{m—n)) 
nature of the triangularization algorithm scales the com¬ 
putational cost drastically with increased model com¬ 
plexity. Furthermore, the need to store a derivative ma¬ 
trix for each iteration results in sizeable memory over¬ 
head when working with large datasets. To mitigate this 
problem, we modified the SRIF algorithm to select a sub¬ 
set of the nominal free parameters during each parameter 
vector adjustment, and to only fit for that subset. We 
tried a variety of subset selection methods, and concluded 
after testing that a “semi-random” mode was the most 
effective. During each parameter vector adjustment, this 
mode randomly selects a fixed number, 6, of parameters 
from the nominal set of parameters {xg} for which the 
condition 

ks<[t^ -\ 
n 

is satisfied, where ks is the total number of times param¬ 
eter Xs has been considered over the course of the entire 
fit, i is the total number of times that the parameter vec¬ 
tor has been adjusted, n is the total number of nominal 
free parameters, and [ J is the round down operator. 

When fittin g for both shape and spin state simultane¬ 
ously (Section l5.3[ ), the spin axis orientation parameters 
were always included in the fit at each parameter vector 
adjustment step. 


A" = 


dR" 

dx 


The algorithm then progresses as described in section 3.1 
with R" replacing R and A" replacing A'. 

5. RESULTS 

We tested our implementation with three different 
types of data. First, we generated simple spherical har¬ 
monic shapes and simulated images with Gaussian noise. 
Second, we used existing shape models of asteroids and 
simulated images with y^-distributed errors, the appro¬ 
priate model for radar noise. Third, we used an actual 
asteroid radar data set. In all cases, we fit the images to 
recover the shapes using SRIF, SPF, and a third-party 
Levenberg-Marquardt algorithm (LM), a standard op¬ 
timizing algorithm which i s used across a wi de variety 


of fields and applications (Press et al. 1992|). Except 
where otherwise noted, our tests did not involve adjust¬ 
ments to parameters controlling the radar scattering law, 
ephemeris corrections, or spin axis orientation. 


5.1. Simulated data with artificial shapes 

These tests consisted of generating an initial basic 
shape (either spherical, oblate ellipsoid, or prolate ellip¬ 
soid), and randomly perturbing the spherical harmonic 
representation of this shape to get a new, non-trivial ob¬ 
ject. 

Simulated range-Doppler images of this object were 
generated, and these images were fit for using the three 
aforementioned algorithms. This test serves as a good 
absolute test of a fitting method, because a solution is 
guaranteed to exist within the framework used (namely, 
a spherical harmonic representation). Figure]^ shows an 
example of the resulting shape when starting with a pro¬ 
late ellipsoid. Three randomly generated objects were 
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created for each of the three basic shapes, for a total of 
nine test cases. 

For each test case, we then generated between 20 and 
30 simulated radar images and added Gaussian noise 
such that pixel values on the target exceed the root- 
mean-square (rms) deviations of the noise by an average 
factor of ^ 5 and a peak factor of ^ 150. These images 
were used to attempt to reconstruct the perturbed shape, 
with the original basic shape (sphere, prolate ellipsoid, 
or oblate ellipsoid) given as the initial condition. This 
process was repeated for each of the nine test cases. 

The three fitting algorithms shared the same starting 
conditions for each test. For each fit, the models com¬ 
prised 121 free parameters (corresponding to the coeffi¬ 
cients of a ten-degree spherical harmonic representation), 
and the simulated images contained a total of 2.4 mil¬ 
lion data points. Stopping criteria were also normalized 
for the three different test types - a fit was considered 
finished if the statistic had not changed to within 
three significant digits after one hour, or twelve hours 
had passed since the fit began, whichever occurred first. 
Time-based stopping criteria - as opposed to iteration- 
based - were chosen in order to account for fundamental 
differences between the algorithms with respect to the 
definition of a single iteration. In addition, fits were al¬ 
lowed to run past the criteria stopping point, and the 
criteria were analyzed and applied afterwards. This was 
to avoid missing a drop-off in in one test type that 
might not appear in another. All times are wall clock 
time. 

The results of these tests (figure ^ indicate that SRIF 
consistently performs better than the currently-used SPF 
algorithm. In addition, SRIF appears to ultimately con¬ 
verge on a lower chi-squared than LM in all cases. 

SRIF also converged on a reduced chi-squared {xled) 
of less than 1.3 (indicating a reasonable approximation 
of the correct model parameters had been found) in eight 
out of the nine tests, while SPF was only able to do so 
in one out of the nine tests. 

5.2. Simulated data with existing asteroid shape models 

We conducted another set of tests using existing shape 
models of asteroids. Three cases were tested - Itokawa, 
the 1999 KW4 primary, and 2000 ET70 (figure |^. As 
opposed to the previous set of tests, these shapes are 
not guaranteed to be well approximated with a spherical 
harmonic representation. However, a best-fit spherical 
harmonic representation can still be found. 

For this set of images we used y^-distributed errors, 
which is the correct noise model for individual images of 
radar echo power. We chose a noise model such that the 
pixel values on the target exceed the rms deviations of 
the noise by an average factor of ^ 2 and a peak factor 
of ^ 60. When multiple images are summed, one can 
often rely on the fact that errors approach normality by 
the central limit theorem, hence the choi ce of Gaussian 
noise for the images tested in section [5T| 

The results of these tests are illustrated in figure 
and an example comparison of a fit synthetic image to 
the simulated image is shown in figure Our implemen¬ 
tation of SRIF clearly performs faster and with higher 
accuracy than both SPF and LM. 

5.3. Spin state 


Jointly solving for spin state and shape is typically 
challenging and time-consuming with the traditional im¬ 
plementation of shape. A common approach is to esti¬ 
mate the spin state as best as possible with rudimentary 
shapes (e.g., ellipsoids or low-order spherical harmonic 
models) in a basic grid search. One can then use the most 
favorable trial values of the spin state to fit a model for 
the physical shape to the observed radar images. Experi¬ 
ence with traditional shape indicates that the algorithm 
rarely deviates much from the initial conditions given on 
the spin state, probably as a result of the one-parameter- 
at-a-time fitting approach. 

Our tests indicate that SRIE is capable of fitting a 
reasonable asteroid shape, even when the initial shape 
and spin state parameters are far away from their opti¬ 
mal values. This advantage likely results from the joint 
estimation of shape and spin parameters. 

Eor example, figure shows the best-fit spherical har¬ 
monics shape, as determined by SRIE, for a set of sim¬ 
ulated images of the asteroid Itokawa. The initial con¬ 
ditions for the shape parameters were a sphere with a 
radius 10% larger than the longest axis of the actual 
shape model. In addition, the initial spin axis was 30 
degrees off from the spin axis with which the data were 
simulated. We repeated these experiments with a variety 
of starting conditions, as well as several different shape 
models, with similar results. 

SRIE’s capacity to fit for both shape and spin state 
parameters can drastically cut down on the total time 
required to obtain an accurate asteroid shape model. 


5.4. Real data: 2000 ET70 

We have run shape with all three fitting algorit hms on 
actual rada r images of the asteroid 2000 ET70 (|Naidu 
et al.| 2013). shape was run initially with an ellipsoid 
model. 4'he starting conditions for this model were such 
that the ellipsoid axes were all equal. The best fit ellip¬ 
soid model \a/h = 1.16, 6/c = 1.13) was then converted 
into a spherical harmonic model with 122 model compo¬ 
nents (corresponding to the coefficients of a ten-degree 
spherical harmonic representation, as well as one over¬ 
all size scaling factor), which was then fit again to the 
data. This process resulted in a final spherical harmonic 
model (figure]^ for the asteroid with a x^ed of 2.1. The 
stopping criterion was a reduction in X^ed than 0.01 
between two consecutive iterations. 

Eor the first stage, SRIE fit a substantially better ellip¬ 
soid than SPE did, although it took about eight minutes 
longer (Table 0. Eor the second stage, SRIE converged 
on a final solution more than two times faster than SPE. 
This further corroborates the results obtained from our 
tests with simulated data. 


6. FUTURE CHANGES 

The addition of SRIE to shape has improved fitting 
performance, but additional changes can still be made to 
allow shape to function optimally with real-world data. 

6.1. Global vs loeal variable partitioning 

The fits discussed in this paper were performed on 
global parameters only - namely, parameters that are 
valid across all data sets associated with the object in 
question. When performing a high-fidelity fit on multi- 
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Figure 2. Example of artifical shape used as a test object to be fitted for. a) The initial shape, a prolate ellipsoid, before any of the 
spherical harmonic parameters have been changed, b) A perturbed prolate ellipsoid. 


Table 1 

Run statistics for SPF and SRIF fits for 2000 ET70 data. 


Model 

Xinitial 

X final 

Runtime (hours) 



SHE 

LM 

SRIF 

SHE 

LM^ 

SRIE 

2000 ET70: Ellipsoid 

3.7 

2.8 

2.7 

2.4 

0.02 

0.84 

0.129 

2000 ET70: Sph. Harm. 

2.4 

2.10 

2.37 

2.10 

2.98 

0.1725 

1.36 


pie data sets, however, it is necessary to take into ac¬ 
count local parameters as well. These are model argu¬ 
ments which apply only at specific points in time. For 
example, while the mean radius and rotation rate of an 
asteroid is a global parameter, the system temperature 
and ephemeris correction parameters on the third day of 
observations are local to the data taken on the third day 
of observations. 

Processing local parameters is less computationally in¬ 
tensive than processing global parameters. The gradients 
of any observables not within a local parameter’s time- 
frame are known a priori to be zero. This greatly re¬ 
duces the number of modelling function calls that must 
be made when considering local parameters. In addi¬ 
tion, the triangularization of a derivative matrix scales 
with the number of non-zero elements. This means that 
while the total number of additional parameters scales 
like the product of the number of datasets with the av¬ 
erage number of local parameters per dataset, the ad¬ 
ditional computation time only scales with the average 
number of local parameters per data set. Because of this, 
adding the capacity for processing local parameters will 
only increase runtime by ^20%. We plan on adding this 
functionality in a future version of shape. 

6.2. Additional fitting methods 

Tests that we have run with the simulated and real 
data have indicated that the y^-space for shape-models 
is not smooth. Figure shows a two-dimensional slice 
of the x^-space for a s^erical harmonic model against 
2000 ET70 data. The multi-valleyed nature of this space 
makes it difficult for local fitting methods to find the 
global minimum. In light of this, global fitting mech¬ 
anisms such as simulated annealing or Markov Chain 
Monte Carlo may be better suited for this problem. 
These methods can be supplemented by a gradient de¬ 
scent method like SRIF. In fact, utilizing a hybrid of 
these methods may prove to be the optimal solution for 
this class of problem. Until such methods are imple¬ 
mented, convergence on a global minimum will be de¬ 
pendent on a good choice of starting conditions. This 
often forces shape modelers to explore a variety of ini¬ 
tial conditions, and identifying such starting conditions 


is not always practical. 

6.3. Additional shape representations 

There are serious drawbacks to using spherical har¬ 
monics to represent the radius of an object at each 
latitude-longitude grid point. Many asteroid shapes are 
poorly approximated by this representation (e.g., the 
1999 KW4 primary) and there are entire classes of shapes 
(e.g., banana) that can not be described at all in this 
fashion. Traditionally, this problem has been solved by 
the use of vertex models, but these shape representations 
typically involve a large number of parameters (i.e., the 
coordinates of three vertices per facet). We are currently 
looking into new representation methods, some of which 
may allow for a greater range of shapes, while at the same 
time cutting down on the number of free parameters. 

7. CONCLUSIONS 

We have added new optimization procedures into the 
asteroid shape modeling software shape, enabling the use 
of a Levenberg-Marquardt algorithm or a Square Root 
Information Filter. We implemented several optimiza¬ 
tions to the SRIF algorithm to increase performance in 
shape inversion problems. Tests on both simulated and 
actual data indicate that our additions allow shape in¬ 
version to proceed more quickly and with better fidelity 
than was previously possible. The SRIF implementation 
also facilitates simultaneous fits of the spin axis orienta¬ 
tion and shape. 















Improved Algorithms for Radar-Based Reconstruction of Asteroid Shapes 


7 


APPENDIX 

NUMERICAL STABILITY 


Issues with numerical stability can arise when multiplying matrices with elements at or near the square root of the 
machine precision. This can lead to erroneous results, or singular matrices for which further calculations (such as the 

matrix inverse) are impossible. _ 

For example ( Bierman, G. J.|[l977 ), consider the case when 


A = 



and 

Then 


W = L 


3 - 2e + 
3-2e 


3-2e 
- 2e + 


(A^WA) = 

Thus, in the case that e is equal to or less than the square root of the machine precision, 

(AWA) = 


2e 3 
2e 3 


2e 

2e 


This matrix is singular, and thus {A'^WA) ^ cannot be computed. This problem is particularly insidious because 
matrix singularity in higher dimensions can be difficult to detect. 
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Figure 3. Results of three fitting algorithms (Sequential Parameter Fit, Levenberg-Marquardt, and Square-Root Information Filter) with 
three artificial shapes (perturbed versions of a sphere, oblate ellipsoid, and prolate ellipsoid). Bold lines indicate fits which converged to a 
^ Dashed lines indicate the assumed future state for fits that had converged on a solution before the 12-hour time frame. 
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25143 Itokawa 2 3 1999 KW4 2 c , 2000 ET70 



1 2 3 4 5 10 15 20 5 10 15 20 

hours hours hours 


Figure 5. Results comparing SRIF, LM, and SPF operating on real asteroid shapes with simulated x^-distributed errors. Dashed lines 
indicate the assumed future state for fits that had converged on a solution. These fits were run without penalty functions. Note that the 
solution arrived at by SPF for 1999 KW4 was a non-physical, pebble-sized asteroid. Avoiding non-physical minima in the x^ space would 
require human intervention to manually tweak the starting conditions. We did not perform such tweaks in order to maintain consistency 
in our tests. 



Figure 6. (a) Example of a simulated input to the shape modeling algorithm. The input is generated by projecting the shape model into 
range-Doppler space at a specific observation epoch and adding random noise, (b) The corresponding synthetic image produced by the 
shape modeling algorithm after fitting for the shape. 
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a) 


b) 

Figure 7. a) The Itokawa shape model that was used to generate simulated radar images (|Ostro et al.|2005|). b) The best-fit SRIF tenth 
degree spherical-harmonics model for those simulated data, using penalty functions. The initial conditions tor the shape parameters were 
a sphere with an offset spin axis latitude and longitude. 
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Figure 8. The best SRIF fit spherical harmonic model for 2000 ET70, which is in good agreement with the model of |Naidu et al.| ( [2013| >. 


^ contour nnap for ET70 



Spherical Harmonic parameter a[0][0] 


1 7.964 
7.319 
6.674 
- 6.029 
-5.384 
- 4.739 
- 4.095 

1 3.450 
2.805 
2.160 


Figure 9. A two-dimensional slice of the space for fitting shape models to radar data. This contour map represents x^ for a spherical 
harmonic model with all parameters fixed except two of the elements of the primary coefficient matrix. Note that the blue regions indicate 
low x ^5 cind that there are several of these regions for a derivative-based optimizer to fall towards, depending on the initial set of starting 
conditions. 
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